home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 2 / AACD 2.iso / AACD / Programming / fpc / inc / lowmath.inc < prev    next >
Text File  |  1998-09-21  |  39KB  |  949 lines
  1. {
  2.     $Id: lowmath.inc,v 1.3 1998/07/01 14:28:32 carl Exp $
  3.     This file is part of the Free Pascal run time library.
  4.     Copyright (c) 1993,97 by Carl-Eric Codere,
  5.     member of the Free Pascal development team
  6.  
  7.     See the file COPYING.FPC, included in this distribution,
  8.     for details about the copyright.
  9.  
  10.     This program is distributed in the hope that it will be useful,
  11.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13.  
  14.  **********************************************************************}
  15. {*************************************************************************}
  16. {  lowmath.inc                                                            }
  17. {  Ported to FPK-Pascal by Carl Eric Codere                               }
  18. {  Terms of use: This source code is freeware.                            }
  19. {*************************************************************************}
  20. {  This inc. implements low-level mathemtical routines  for the motorola  }
  21. {  68000 family of processors.                                            }
  22. {*************************************************************************}
  23. { single floating point routines taken from GCC 2.5.2 Atari compiler      }
  24. { library source.                                                         }
  25. {  Original credits:                                                      }
  26. {    written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).                   }
  27. {      Based on a 80x86 floating point packet from comp.os.minix,         }
  28. {        written by P.Housel                                              }
  29. {    Patched by Olaf Flebbe (flebbe@tat.physik.uni-tuebingen.de)          }
  30. {    Revision by michal 05-93 (ntomczak@vm.ucs.ualberta.ca)               }
  31. {*************************************************************************}
  32. {--------------------------------------------------------------------}
  33. { LEFT TO DO:                                                        }
  34. { o Add support for FPU if present.                                  }
  35. { o Verify if single comparison works in all cases.                  }
  36. { o Add support for  NANs in SINGLE_CMP                              }
  37. { o Add comp (80-bit) multiplication,addition,substract,division,    }
  38. {    shift.                                                          }
  39. { o Add stack checking for the routines which use the stack.         }
  40. {    (This will probably have to be done in the code generator).     }
  41. {--------------------------------------------------------------------}
  42.  
  43.  
  44.  
  45. Procedure Single_Norm;Assembler;
  46. {--------------------------------------------}
  47. { Low-level routine to normalize single  e   }
  48. {  IEEE floating point values. Never called  }
  49. {  directly.                                 }
  50. { On Exit:                                   }
  51. {      d0 = result.                          }
  52. {  Registers destroyed: d0,d1                }
  53. {--------------------------------------------}
  54. Asm
  55.    tst.l    d4             { rounding and u.mant == 0 ?    }
  56.    bne      @normlab1
  57.    tst.b    d1
  58.    beq      @retzok
  59. @normlab1:
  60.    clr.b    d2             { "sticky byte"                 }
  61. @normlab3:
  62.    move.l   #$ff000000,d5
  63. @normlab4:
  64.    tst.w    d0             { divide (shift)                }
  65.    ble      @normlab2      {  denormalized number          }
  66.    move.l   d4,d3
  67.    and.l    d5,d3          {  or until no bits above 23    }
  68.    beq      @normlab5
  69. @normlab2:
  70.    addq.w   #1,d0          { increment exponent            }
  71.    lsr.l    #1,d4
  72.    or.b     d1,d2          { set "sticky"                  }
  73.    roxr.b   #1,d1          { shift into rounding bits      }
  74.    bra      @normlab4
  75. @normlab5:
  76.    and.b    #1,d2
  77.    or.b     d2,d1          { make least sig bit "sticky"   }
  78.    asr.l    #1,d5          { #0xff800000 -> d5             }
  79. @normlab6:
  80.    move.l   d4,d3          { multiply (shift) until        }
  81.    and.l    d5,d3          { one in "implied" position     }
  82.    bne      @normlab7
  83.    subq.w   #1,d0          { decrement exponent            }
  84.    beq      @normlab7      {  too small. store as denormalized number }
  85.    add.b    d1,d1          { some doubt about this one *              }
  86.    addx.l   d4,d4
  87.    bra      @normlab6
  88. @normlab7:
  89.    tst.b    d1             { check rounding bits          }
  90.    bge      @normlab9      { round down - no action neccessary }
  91.    neg.b    d1
  92.    bvc      @normlab8      { round up                     }
  93.    move.w   d4,d1          { tie case - round to even     }
  94.                            {   dont need rounding bits any more }
  95.    and.w    #1,d1          { check if even                }
  96.    beq      @normlab9      {   mantissa is even - no action necessary  }
  97.                            { fall through                 }
  98. @normlab8:
  99.    clr.w    d1             { zero rounding bits            }
  100.    add.l    #1,d4
  101.    tst.w    d0
  102.    bne      @normlab10     { renormalize if number was denormalized   }
  103.    add.w    #1,d0          { correct exponent for denormalized numbers }
  104.    bra      @normlab3
  105. @normlab10:
  106.    move.l   d4,d3          { check for rounding overflow              }
  107.    asl.l    #1,d5          { #0xff000000 -> d5                        }
  108.    and.l    d5,d3
  109.    bne      @normlab4      { go back and renormalize                  }
  110. @normlab9:
  111.    tst.l    d4             { check if normalization caused an underflow }
  112.    beq      @retz
  113.    tst.w    d0             { check for exponent overflow or underflow  }
  114.    blt      @retz
  115.    cmp.w    #255,d0
  116.    bge      @oflow
  117.  
  118.    lsl.w    #8,d0          { re-position exponent - one bit too high }
  119.    lsl.w    #1,d2          { get X bit                               }
  120.    roxr.w   #1,d0          { shift it into sign position             }
  121.    swap     d0             { map to upper word                       }
  122.    clr.w    d0
  123.    and.l    #$7fffff,d4    { top mantissa bits                       }
  124.    or.l     d4,d0          { insert exponent and sign                }
  125.    movem.l  (sp)+,d2-d5
  126.    rts
  127.  
  128. @retz:
  129.    { handling underflow should be done here... }
  130.    { by default simply return 0 as retzok...   }
  131. @retzok:
  132.    moveq.l   #0,d0
  133.    lsl.w     #1,d2
  134.    roxr.l    #1,d0          { sign of 0 is the same as of d2          }
  135.    movem.l   (sp)+,d2-d5
  136.    rts
  137.  
  138. @oflow:
  139.    move.l    #$7f800000,d0  { +infinity as proposed by IEEE         }
  140.  
  141.    tst.w     d2             { transfer sign                         }
  142.    bge       @ofl_clear     { (mjr++)                               }
  143.    bset      #31,d0         {                                       }
  144. @ofl_clear:
  145.    or.b      #2,ccr         { set overflow flag.                    }
  146.    movem.l   (sp)+,d2-d5
  147.    rts
  148. end;
  149.  
  150.  
  151. Procedure Single_AddSub; Assembler;
  152. {--------------------------------------------}
  153. { Low-level routine to add/subtract single   }
  154. {  IEEE floating point values. Never called  }
  155. {  directly.                                 }
  156. { On Exit:                                   }
  157. {      d0 = result -- from normalize routine }
  158. {      Flags : V set if overflow.            }
  159. {              on underflow d0 = 0           }
  160. {  Registers destroyed: d0,d1                }
  161. {--------------------------------------------}
  162. Asm
  163. {--------------------------------------------}
  164. { On Entry:                                  }
  165. {      d1-d0 = single values to subtract.    }
  166. {--------------------------------------------}
  167. XDEF SINGLE_SUB
  168.    eor.l      #$80000000,d0         {   reverse sign of v }
  169. {--------------------------------------------}
  170. { On Entry:                                  }
  171. {      d0, d1 = single values to add.        }
  172. {--------------------------------------------}
  173. XDEF SINGLE_ADD
  174.    movem.l d2-d5,-(sp)              { save registers      }
  175.    move.l   d0,d4                   { d4 = d0 = v         }
  176.    move.l   d1,d5                   { d5 = d1 = u         }
  177.  
  178.    move.l  #$7fffff,d3
  179.    move.l  d5,d0                    { d0 = u.exp          }
  180.    move.l  d5,d2                    { d2.h = u.sign       }
  181.    swap       d0
  182.    move.w  d0,d2                    { d2 = u.sign         }
  183.    and.l   d3,d5                    { remove exponent from u.mantissa }
  184.  
  185.    move.l  d4,d1                    { d1 = v.exp          }
  186.    and.l   d3,d4                    { remove exponent from v.mantissa }
  187.    swap    d1
  188.    eor.w   d1,d2                    { d2 = u.sign ^ v.sign (in bit 15)}
  189.    clr.b   d2                       { we will use the lowest byte as a flag }
  190.    moveq.l #15,d3
  191.    bclr    d3,d1                    { kill sign bit u.exp             }
  192.    bclr    d3,d0                    { kill sign bit u.exp             }
  193.    btst    d3,d2                    { same sign for u and v?          }
  194.    beq     @slabel1
  195.    cmp.l   d0,d1                    { different signs - maybe x - x ? }
  196.    seq     d2                       { set 'cancellation' flag         }
  197. @slabel1:
  198.    lsr.w   #7,d0                    { keep here exponents only        }
  199.    lsr.w   #7,d1
  200. {--------------------------------------------------------------------}
  201. {          Now perform testing of NaN and infinities                 }
  202. {--------------------------------------------------------------------}
  203.    moveq.l #-1,d3
  204.    cmp.b   d3,d0
  205.    beq     @alabel1
  206.    cmp.b   d3,d1
  207.    bne     @nospec
  208.    bra     @alabel2
  209. {--------------------------------------------------------------------}
  210. {                        u is special.                               }
  211. {--------------------------------------------------------------------}
  212. @alabel1:
  213.    tst.b      d2
  214.    bne        @retnan      { cancellation of specials -> NaN }
  215.    tst.l      d5
  216.    bne        @retnan      { arith with Nan gives always NaN }
  217.  
  218.    addq.w     #4,a0        { here is an infinity             }
  219.    cmp.b      d3,d1
  220.    bne        @alabel3     { skip check for NaN if v not special }
  221. {--------------------------------------------------------------------}
  222. {                        v is special.                               }
  223. {--------------------------------------------------------------------}
  224. @alabel2:
  225.    tst.l           d4
  226.    bne        @retnan
  227. @alabel3:
  228.    move.l     (a0),d0
  229.    bra        @return
  230. {--------------------------------------------------------------------}
  231. {                     Return a quiet nan                             }
  232. {--------------------------------------------------------------------}
  233. @retnan:
  234.    moveq.l    #-1,d0
  235.    lsr.l      #1,d0                { 0x7fffffff -> d0        }
  236.    bra        @return
  237. { Ok, no inifinty or NaN involved.. }
  238. @nospec:
  239.    tst.b           d2
  240.    beq        @alabel4
  241.    moveq.l    #0,d0                { x - x hence we always return +0 }
  242. @return:
  243.    movem.l    (sp)+,d2-d5
  244.    rts
  245.  
  246. @alabel4:
  247.    moveq.l    #23,d3
  248.    bset       d3,d5           { restore implied leading "1"   }
  249.    tst.w      d0              { check for zero exponent - no leading "1" }
  250.    bne        @alabel5
  251.    bclr       d3,d5           { remove it                     }
  252.    addq.w     #1,d0           { "normalize" exponent          }
  253. @alabel5:
  254.    bset       d3,d4           { restore implied leading "1"   }
  255.    tst.w      d1              { check for zero exponent - no leading "1"  }
  256.    bne        @alabel6
  257.    bclr       d3,d4           { remove it                     }
  258.    addq.w  #1,d1              { "normalize" exponent          }
  259. @alabel6:
  260.    moveq.l    #0,d3           { (put initial zero rounding bits in d3) }
  261.    neg.w      d1              { d1 = u.exp - v.exp            }
  262.    add.w      d0,d1
  263.    beq        @alabel8      { exponents are equal - no shifting neccessary }
  264.    bgt        @alabel7      { not equal but no exchange neccessary         }
  265.    exg        d4,d5                { exchange u and v }
  266.    sub.w      d1,d0                { d0 = u.exp - (u.exp - v.exp) = v.exp }
  267.    neg.w      d1
  268.    tst.w      d2              { d2.h = u.sign ^ (u.sign ^ v.sign) = v.sign }
  269.    bpl        @alabel7
  270.    bchg       #31,d2
  271. @alabel7:
  272.    cmp.w      #26,d1       { is u so much bigger that v is not }
  273.    bge        @alabel9      { significant ?                     }
  274. {--------------------------------------------------------------------}
  275. { shift mantissa left two digits, to allow cancellation of           }
  276. { most significant digit, while gaining an additional digit for      }
  277. { rounding.                                                          }
  278. {--------------------------------------------------------------------}
  279.    moveq.l #1,d3
  280. @alabel10:
  281.    add.l           d5,d5
  282.    subq.w  #1,d0              { decrement exponent            }
  283.    subq.w  #1,d1              { done shifting altogether ?    }
  284.    dbeq    d3,@alabel10        { loop if still can shift u.mant more }
  285.    moveq.l    #0,d3
  286.  
  287.    cmp.w      #16,d1          { see if fast rotate possible    }
  288.    blt        @alabel11
  289.    or.w       d4,d3           { set rounding bits              }
  290.    clr.w      d4
  291.    swap       d4
  292.    subq.w  #8,d1
  293.    subq.w  #8,d1
  294.    bra      @alabel11
  295.  
  296. @alabel12:
  297.    move.b    d4,d2
  298.    and.b      #1,d2
  299.    or.b       d2,d3
  300.    lsr.l      #1,d4               { shift v.mant right the rest of the way }
  301. @alabel11:
  302.    dbra    d1,@alabel12           { loop                                   }
  303.  
  304. @alabel8:
  305.    tst.w      d2                  { are the signs equal ?                  }
  306.    bpl        @alabel13           { yes, no negate necessary               }
  307.  
  308.  
  309.    tst.w      d3                  { negate rounding bits and v.mant        }
  310.    beq        @alabel14
  311.    addq.l  #1,d4
  312. @alabel14:
  313.    neg.l           d4
  314.  
  315. @alabel13:
  316.    add.l      d4,d5                { u.mant = u.mant + v.mant              }
  317.    bcs        @alabel9             { needn not negate                      }
  318.    tst.w      d2                   { opposite signs ?                      }
  319.    bpl        @alabel9             { do not need to negate result          }
  320.  
  321.    neg.l      d5
  322.    not.l      d2                   { switch sign                           }
  323. @alabel9:
  324.    move.l  d5,d4                   { move result for normalization         }
  325.    clr.l      d1
  326.    tst.l      d3
  327.    beq        @alabel15
  328.    moveq.l #-1,d1
  329. @alabel15:
  330.    swap    d2                      { put sign into d2 (exponent is in d0)  }
  331.    jmp     SINGLE_NORM             { leave registers on stack for norm_sf  }
  332. end;
  333.  
  334.  
  335. Procedure Single_Mul;Assembler;
  336. {--------------------------------------------}
  337. { Low-level routine to multiply two single   }
  338. {  IEEE floating point values. Never called  }
  339. {  directly.                                 }
  340. { Om Entry:                                  }
  341. {      d0,d1 = values to multiply            }
  342. { On Exit:                                   }
  343. {      d0 = result.                          }
  344. {  Registers destroyed: d0,d1                }
  345. { stack space used (and restored): 8 bytes.  }
  346. {--------------------------------------------}
  347. Asm
  348. XDEF SINGLE_MUL
  349.    movem.l  d2-d5,-(sp)
  350.    move.l   d0,d4       { d4 = v                   }
  351.    move.l   d1,d5       { d5 = u                   }
  352.  
  353.    move.l   #$7fffff,d3
  354.    move.l   d5,d0       { d0 = u.exp               }
  355.    and.l    d3,d5       { remove exponent from u.mantissa }
  356.    swap     d0
  357.    move.w   d0,d2       { d2 = u.sign              }
  358.  
  359.    move.l   d4,d1       { d1 = v.exp               }
  360.    and.l    d3,d4       { remove exponent from v.mantissa }
  361.    swap     d1
  362.    eor.w    d1,d2       { d2 = u.sign ^ v.sign (in bit 15)}
  363.  
  364.    moveq.l  #15,d3
  365.    bclr     d3,d0       { kill sign bit                   }
  366.    bclr     d3,d1       { kill sign bit                   }
  367.    tst.l    d0          { test if one of factors is 0     }
  368.    beq      @mlabel1
  369.    tst.l    d1
  370. @mlabel1:
  371.    seq      d2         { 'one of factors is 0' flag in the lowest byte }
  372.    lsr.w    #7,d0      { keep here exponents only         }
  373.    lsr.w    #7,d1
  374.  
  375. {--------------------------------------------------------------------}
  376. {          Now perform testing of NaN and infinities                 }
  377. {--------------------------------------------------------------------}
  378.    moveq.l  #-1,d3
  379.    cmp.b    d3,d0
  380.    beq      @mlabel2
  381.    cmp.b    d3,d1
  382.    bne      @mnospec
  383.    bra      @mlabel3
  384. {--------------------------------------------------------------------}
  385. {                   first operand is special                         }
  386. {--------------------------------------------------------------------}
  387. @mlabel2:
  388.    tst.l    d5         { is it NaN?                       }
  389.    bne      @mretnan
  390. @mlabel3:
  391.    tst.b    d2         { 0 times special or special times 0 ? }
  392.    bne      @mretnan   { yes -> NaN                       }
  393.    cmp.b    d3,d1      { is the other special ?           }
  394.    beq      @mlabel4   { maybe it is NaN                  }
  395. {--------------------------------------------------------------------}
  396. {                   Return infiny with correct sign                  }
  397. {--------------------------------------------------------------------}
  398. @mretinf:
  399.    move.l   #$ff000000,d0  { we will return #0xff800000 or #0x7f800000 }
  400.    lsl.w    #1,d2
  401.    roxr.l   #1,d0      { shift in high bit as given by d2 }
  402. @mreturn:
  403.    movem.l  (sp)+,d2-d5
  404.    rts
  405.  
  406. {--------------------------------------------------------------------}
  407. {                        v is special.                               }
  408. {--------------------------------------------------------------------}
  409. @mlabel4:
  410.    tst.l    d4          { is this NaN?                    }
  411.    beq      @mretinf    { we know that the other is not zero }
  412. @mretnan:
  413.    moveq.l  #-1,d0
  414.    lsr.l    #1,d0       { 0x7fffffff -> d0                }
  415.    bra      @mreturn
  416. {--------------------------------------------------------------------}
  417. {                    End of NaN and Inf                              }
  418. {--------------------------------------------------------------------}
  419. @mnospec:
  420.    tst.b    d2          { not needed - but we can waste two instr. }
  421.    bne      @mretzz     { return signed 0 if one of factors is 0   }
  422.    moveq.l  #23,d3
  423.    bset     d3,d5       { restore implied leading "1"              }
  424.    subq.w   #8,sp       { multiplication accumulator               }
  425.    tst.w    d0          { check for zero exponent - no leading "1" }
  426.    bne      @mlabel5
  427.    bclr     d3,d5       { remove it                                }
  428.    addq.w   #1,d0       { "normalize" exponent                     }
  429. @mlabel5:
  430.    tst.l   d5
  431.    beq     @mretz       { multiplying zero                         }
  432.  
  433.    moveq.l #23,d3
  434.    bset    d3,d4        { restore implied leading "1"              }
  435.    tst.w   d1           { check for zero exponent - no leading "1" }
  436.    bne     @mlabel6
  437.    bclr    d3,d4        { remove it                                }
  438.    addq.w  #1,d1        { "normalize" exponent                     }
  439. @mlabel6:
  440.    tst.l   d4
  441.    beq     @mretz       { multiply by zero                         }
  442.  
  443.    add.w   d1,d0        { add exponents,                           }
  444.    sub.w   #BIAS4+16-8,d0 { remove excess bias, acnt for repositioning }
  445.  
  446.    clr.l   (sp)         { initialize 64-bit product to zero        }
  447.    clr.l   4(sp)
  448. {--------------------------------------------------------------------}
  449. { see Knuth, Seminumerical Algorithms, section 4.3. algorithm M      }
  450. {--------------------------------------------------------------------}
  451.    move.w  d4,d3
  452.    mulu.w  d5,d3        { mulitply with bigit from multiplier      }
  453.    move.l  d3,4(sp)     { store into result                        }
  454.  
  455.    move.l  d4,d3
  456.    swap    d3
  457.    mulu.w  d5,d3
  458.    add.l   d3,2(sp)     { add to result                            }
  459.  
  460.    swap    d5           { [TOP 8 BITS SHOULD BE ZERO !]            }
  461.  
  462.    move.w  d4,d3
  463.    mulu.w  d5,d3        { mulitply with bigit from multiplier      }
  464.    add.l   d3,2(sp)     { store into result (no carry can occur here) }
  465.  
  466.    move.l  d4,d3
  467.    swap    d3
  468.    mulu.w  d5,d3
  469.    add.l   d3,(sp)      { add to result                            }
  470.             { [TOP 16 BITS SHOULD BE ZERO !] }
  471.    movem.l 2(sp),d4-d5  { get the 48 valid mantissa bits           }
  472.    clr.w   d5           { (pad to 64)                              }
  473.  
  474.    move.l  #$0000ffff,d3
  475. @mlabel7:
  476.    cmp.l   d3,d4        { multiply (shift) until                  }
  477.    bhi     @mlabel8     {  1 in upper 16 result bits              }
  478.    cmp.w   #9,d0        { give up for denormalized numbers        }
  479.    ble     @mlabel8
  480.    swap    d4         { (we''re getting here only when multiplying }
  481.    swap    d5         {  with a denormalized number; there''s an   }
  482.    move.w  d5,d4      {  eventual loss of 4 bits in the rounding   }
  483.    clr.w   d5         {  byte -- what a pity 8-)                   }
  484.    subq.w  #8,d0      { decrement exponent                         }
  485.    subq.w  #8,d0
  486.    bra     @mlabel7
  487. @mlabel8:
  488.    move.l  d5,d1      { get rounding bits                          }
  489.    rol.l   #8,d1
  490.    move.l  d1,d3      { see if sticky bit should be set            }
  491.    and.l   #$ffffff00,d3
  492.    beq     @mlabel9
  493.    or.b    #1,d1      { set "sticky bit" if any low-order set      }
  494. @mlabel9:
  495.    addq.w  #8,sp      { remove accumulator from stack              }
  496.    jmp     SINGLE_NORM{ (result in d4)                             }
  497.  
  498. @mretz:
  499.    addq.w  #8,sp      { release accumulator space                  }
  500. @mretzz:
  501.    moveq.l #0,d0      { save zero as result                        }
  502.    lsl.w   #1,d2      { and set it sign as for d2                  }
  503.    roxr.l  #1,d0
  504.    movem.l (sp)+,d2-d5
  505.    rts                { no normalizing neccessary                  }
  506. end;
  507.  
  508.  
  509. Procedure Single_Div;Assembler;
  510. {--------------------------------------------}
  511. { Low-level routine to dividr two single     }
  512. {  IEEE floating point values. Never called  }
  513. {  directly.                                 }
  514. { Om Entry:                                  }
  515. {      d1/d0 = u/v = operation to perform.   }
  516. { On Exit:                                   }
  517. {      d0 = result.                          }
  518. {  Registers destroyed: d0,d1                }
  519. { stack space used (and restored): 8 bytes.  }
  520. {--------------------------------------------}
  521. ASM
  522. XDEF SINGLE_DIV
  523.    { u = d1 = dividend }
  524.    { v = d0 = divisor  }
  525.    tst.l  d0              { check if divisor is 0                  }
  526.    bne    @dno_exception
  527.  
  528.    move.l #$7f800000,d0
  529.    btst   #31,d1          { transfer sign of dividend              }
  530.    beq    @dclear
  531.    bset   #31,d0
  532. @dclear:
  533.    rts
  534.  
  535. @dno_exception:
  536.    move.l  d1,d4          { d4 = u, d5 = v                         }
  537.    move.l  d0,d5
  538.    movem.l d2-d5,-(sp)    { save registers                         }
  539.  
  540.    move.l  #$7fffff,d3
  541.    move.l  d4,d0          { d0 = u.exp                             }
  542.    and.l   d3,d4          { remove exponent from u.mantissa        }
  543.    swap    d0
  544.    move.w  d0,d2          { d2 = u.sign                            }
  545.  
  546.    move.l  d5,d1          { d1 = v.exp                             }
  547.    and.l   d3,d5          { remove exponent from v.mantissa        }
  548.    swap    d1
  549.    eor.w   d1,d2          { d2 = u.sign ^ v.sign (in bit 15)       }
  550.  
  551.    moveq.l #15,d3
  552.    bclr    d3,d0          { kill sign bit                          }
  553.    bclr    d3,d1          { kill sign bit                          }
  554.    lsr.w   #7,d0
  555.    lsr.w   #7,d1
  556.  
  557.    moveq.l #-1,d3
  558.    cmp.b   d3,d0          { comparison with #0xff                  }
  559.    beq     @dlabel1       { u == NaN ;; u== Inf                    }
  560.    cmp.b   d3,d1
  561.    beq     @dlabel2       { v == NaN ;; v == Inf                   }
  562.    tst.b   d0
  563.    bne     @dlabel4       { u not zero nor denorm                  }
  564.    tst.l   d4
  565.    beq     @dlabel3       { 0/ ?                                   }
  566.  
  567. @dlabel4:
  568.    tst.w   d1
  569.    bne     @dnospec
  570.  
  571.    tst.l   d5
  572.    bne     @dnospec
  573.    bra     @dretinf       { x/0 -> +/- Inf                         }
  574.  
  575. @dlabel1:
  576.    tst.l   d4             { u == NaN ?                             }
  577.    bne     @dretnan       { NaN/ x                                 }
  578.    cmp.b   d3,d1
  579.    beq     @dretnan       { Inf/Inf or Inf/NaN                     }
  580. {  bra     dretinf      ; Inf/x ; x != Inf && x != NaN             }
  581. {--------------------------------------------------------------------}
  582. {                  Return infinity with correct sign.                }
  583. {--------------------------------------------------------------------}
  584. @dretinf:
  585.    move.l  #$ff000000,d0
  586.    lsl.w   #1,d2
  587.    roxr.l  #1,d0          { shift in high bit as given by d2       }
  588. @dreturn:
  589.    movem.l (sp)+,d2-d5
  590.    rts
  591.  
  592. @dlabel2:
  593.    tst.l   d5
  594.    bne     @dretnan       { x/NaN                                  }
  595. {   bra    dretzero   ; x/Inf -> +/- 0                             }
  596. {--------------------------------------------------------------------}
  597. {                  Return correct signed zero.                       }
  598. {--------------------------------------------------------------------}
  599. @dretzero:
  600.    moveq.l #0,d0          { zero destination                       }
  601.    lsl.w   #1,d2          { set X bit accordingly                  }
  602.    roxr.l  #1,d0
  603.    bra     @dreturn
  604.  
  605. @dlabel3:
  606.    tst.w   d1
  607.    bne     @dretzero      { 0/x ->+/- 0                            }
  608.    tst.l   d4
  609.    bne     @dretzero      { 0/x                                    }
  610. {   bra    dretnan         0/0                                     }
  611. {--------------------------------------------------------------------}
  612. {                       Return NotANumber                            }
  613. {--------------------------------------------------------------------}
  614. @dretnan:
  615.    move.l  d3,d0          { d3 contains 0xffffffff                 }
  616.    lsr.l   #1,d0
  617.    bra     @dreturn
  618. {--------------------------------------------------------------------}
  619. {                      End of Special Handling                       }
  620. {--------------------------------------------------------------------}
  621. @dnospec:
  622.    moveq.l #23,d3
  623.    bset    d3,d4          { restore implied leading "1"            }
  624.    tst.w   d0             { check for zero exponent - no leading "1" }
  625.    bne     @dlabel5
  626.    bclr    d3,d4          { remove it                              }
  627.    add.w   #1,d0          { "normalize" exponent                   }
  628. @dlabel5:
  629.    tst.l   d4
  630.    beq     @dretzero      { dividing zero                          }
  631.  
  632.    bset    d3,d5          { restore implied leading "1"            }
  633.    tst.w   d1             { check for zero exponent - no leading "1"}
  634.    bne     @dlabel6
  635.    bclr    d3,d5          { remove it                              }
  636.    add.w   #1,d1          { "normalize" exponent                   }
  637. @dlabel6:
  638.  
  639.    sub.w   d1,d0          { subtract exponents,                    }
  640.    add.w   #BIAS4-8+1,d0  { add bias back in, account for shift    }
  641.    add.w   #34,d0     { add loop offset, +2 for extra rounding bits}
  642.                       { for denormalized numbers (2 implied by dbra)}
  643.    move.w  #27,d1     { bit number for "implied" pos (+4 for rounding)}
  644.    moveq.l #-1,d3     { zero quotient (for speed a one''s complement) }
  645.    sub.l   d5,d4      { initial subtraction, u = u - v                }
  646. @dlabel7:
  647.    btst    d1,d3      { divide until 1 in implied position            }
  648.    beq     @dlabel9
  649.  
  650.    add.l   d4,d4
  651.    bcs     @dlabel8   { if carry is set, add, else subtract           }
  652.  
  653.    addx.l  d3,d3      { shift quotient and set bit zero               }
  654.    sub.l   d5,d4      { subtract, u = u - v                           }
  655.    dbra    d0,@dlabel7      { give up if result is denormalized        }
  656.    bra     @dlabel9
  657. @dlabel8:
  658.    addx.l  d3,d3      { shift quotient and clear bit zero             }
  659.    add.l   d5,d4      { add (restore), u = u + v                      }
  660.    dbra    d0,@dlabel7      { give up if result is denormalized        }
  661. @dlabel9:
  662.    subq.w  #2,d0      { remove rounding offset for denormalized nums  }
  663.    not.l   d3         { invert quotient to get it right               }
  664.  
  665.    clr.l   d1         { zero rounding bits                            }
  666.    tst.l   d4         { check for exact result                        }
  667.    beq     @dlabel10
  668.    moveq.l #-1,d1     { prevent tie case                              }
  669. @dlabel10:
  670.    move.l  d3,d4      { save quotient mantissa                        }
  671.    jmp     SINGLE_NORM{ (registers on stack removed by norm_sf)       }
  672. end;
  673.  
  674.  
  675. Procedure Single_Cmp; Assembler;
  676. {--------------------------------------------}
  677. { Low-level routine to compare single two    }
  678. {  single point values..                     }
  679. {  Never called directly.                    }
  680. { On Entry:                                  }
  681. {      d1 and d0 Values to compare           }
  682. {      d0 = first operand                    }
  683. { On Exit:                                   }
  684. {      Flags according to result             }
  685. {  Registers destroyed: d0,d1                }
  686. {--------------------------------------------}
  687. Asm
  688. XDEF  SINGLE_CMP
  689.    tst.l     d0               { check sign bit             }
  690.    bpl       @cmplab1
  691.    neg.l     d0               { negate                     }
  692.    bchg      #31,d0           { toggle sign bit            }
  693. @cmplab1:
  694.    tst.l     d1               { check sign bit             }
  695.    bpl       @cmplab2
  696.    neg.l     d1               { negate                     }
  697.    bchg      #31,d1           { toggle sign bit            }
  698. @cmplab2:
  699.    cmp.l     d0,d1            { compare...                 }
  700.    rts
  701. end;
  702.  
  703.  
  704.  
  705. Procedure LongMul;Assembler;
  706. {--------------------------------------------}
  707. { Low-level routine to multiply two signed   }
  708. {  32-bit values. Never called directly.     }
  709. { On entry: d1,d0 = 32-bit signed values to  }
  710. {           multiply.                        }
  711. { On Exit:                                   }
  712. {      d0 = result.                          }
  713. {  Registers destroyed: d0,d1                }
  714. {  stack space used and restored: 10 bytes   }
  715. {--------------------------------------------}
  716. Asm
  717. XDEF LONGMUL
  718.              cmp.b   #2,Test68000    { Are we on a 68020+ cpu                  }
  719.              blt     @Lmulcontinue
  720.              muls.l  d1,d0           { yes, then directly mul...               }
  721.              rts                     { return... result in d0                  }
  722. @Lmulcontinue:
  723.              move.l    d2,a0         { save registers                          }
  724.              move.l    d3,a1
  725.  
  726.              move.l    d0,-(sp)
  727.              move.l    d1,-(sp)
  728.  
  729.              movem.w   (sp)+,d0-d3   { u = d0-d1, v = d2-d3                    }
  730.  
  731.              move.w    d0,-(sp)      { sign flag                               }
  732.              bpl       @LMul1        { is u negative ?                         }
  733.              neg.w     d1            { yes, force it positive                  }
  734.              negx.w    d0
  735. @LMul1:
  736.              tst.w     d2            { is v negative ?                         }
  737.              bpl       @LMul2
  738.              neg.w     d3            { yes, force it positive ...              }
  739.              negx.w    d2
  740.              not.w     (sp)          {  ... and modify flag word               }
  741. @LMul2:
  742.              ext.l     d0            { u.h <> 0 ?                              }
  743.              beq       @LMul3
  744.              mulu.w    d3,d0         { r  = v.l * u.h                          }
  745. @LMul3:
  746.              tst.w     d2            { v.h <> 0 ?                              }
  747.              beq       @LMul4
  748.              mulu.w    d1,d2         { r += v.h * u.l                          }
  749.              add.w     d2,d0
  750. @LMul4:
  751.              swap      d0
  752.              clr.w     d0
  753.              mulu.w    d3,d1        { r += v.l * u.l                           }
  754.              add.l     d1,d0
  755.              move.l    a1,d3
  756.              move.l    a0,d2
  757.              tst.w     (sp)+        { should the result be negated ?           }
  758.              bpl       @LMul5       { no, just return                          }
  759.              neg.l     d0           { else r = -r                              }
  760. @LMul5:
  761.              rts
  762. end;
  763.  
  764.  
  765.  
  766. Procedure Long2Single;Assembler;
  767. {--------------------------------------------}
  768. { Low-level routine to convert a longint     }
  769. {  to a single floating point value.         }
  770. { On entry: d0 = longint value to convert.   }
  771. { On Exit:                                   }
  772. {      d0 = single IEEE value                }
  773. {  Registers destroyed: d0,d1                }
  774. {  stack space used and restored:  8 bytes   }
  775. {--------------------------------------------}
  776. Asm
  777. XDEF LONG2SINGLE
  778.    movem.l  d2-d5,-(sp)  { save registers to make norm_sf happy}
  779.  
  780.    move.l   d0,d4        { prepare result mantissa          }
  781.    move.w   #BIAS4+32-8,d0   { radix point after 32 bits    }
  782.    move.l   d4,d2        { set sign flag                    }
  783.    bge      @l2slabel1   { nonnegative                      }
  784.    neg.l    d4           { take absolute value              }
  785. @l2slabel1:
  786.    swap     d2           { follow SINGLE_NORM conventions   }
  787.    clr.w    d1           { set rounding = 0                 }
  788.    jmp      SINGLE_NORM
  789. end;
  790.  
  791.  
  792. Procedure LongDiv; Assembler;
  793. {--------------------------------------------}
  794. { Low-level routine to do signed long        }
  795. { division.                                  }
  796. { On entry: d0/d1 operation to perform       }
  797. { On Exit:                                   }
  798. {      d0 = quotient                         }
  799. {      d1 = remainder                        }
  800. {  Registers destroyed: d0,d1,d6             }
  801. {  stack space used and restored: 10 bytes   }
  802. {--------------------------------------------}
  803. asm
  804. XDEF LONGDIV
  805.    cmp.b       #2,Test68000  { can we use divs ?                  }
  806.    blt         @continue
  807.    tst.l       d1
  808.    beq         @zerodiv2
  809.    move.l      d1,d6
  810.    clr.l       d1           { clr                                 }
  811.    tst.l       d0           { check sign of d0                    }
  812.    bpl         @posdiv
  813.    move.l      #$ffffffff,d1{ sign extend into  d1                }
  814. @posdiv:
  815.    divsl.l     d6,d1:d0
  816.    rts
  817. @continue:
  818.  
  819.    move.l      d2,a0      { save registers                        }
  820.    move.l      d3,a1
  821.    move.l      d4,-(sp)   { divisor = d1 = d4                     }
  822.    move.l      d5,-(sp)   { divident = d0 = d5                    }
  823.  
  824.    move.l      d1,d4      { save divisor                          }
  825.    move.l      d0,d5      { save dividend                         }
  826.  
  827.    clr.w       -(sp)      { sign flag                             }
  828.  
  829.    clr.l       d0         { prepare result                        }
  830.    move.l      d4,d2      { get divisor                           }
  831.    beq         @zerodiv   { divisor = 0 ?                         }
  832.    bpl         @LDiv1     { divisor < 0 ?                         }
  833.    neg.l       d2         { negate it                             }
  834.    not.w       (sp)       { remember sign                         }
  835. @LDiv1:
  836.    move.l      d5,d1      { get dividend                          }
  837.    bpl         @LDiv2     { dividend < 0 ?                        }
  838.    neg.l       d1         { negate it                             }
  839.    not.w       (sp)       { remember sign                         }
  840. @LDiv2:
  841. {;== case 1) divident < divisor}
  842.    cmp.l       d2,d1      { is divident smaller then divisor ?    }
  843.    bcs         @LDiv7     { yes, return immediately               }
  844. {;== case 2) divisor has <= 16 significant bits}
  845.    move.l      d4,d6      { put divisor in d6 register            }
  846.    lsr.l       #8,d6      { rotate into low word                  }
  847.    lsr.l       #8,d6
  848.    tst.l       d6
  849.    bne         @LDiv3     { divisor has only 16 bits              }
  850.    move.w      d1,d3      { save dividend                         }
  851.    clr.w       d1         { divide dvd.h by dvs                   }
  852.    swap        d1
  853.    beq         @LDiv4     { (no division necessary if dividend zero)}
  854.    divu        d2,d1
  855. @LDiv4:
  856.    move.w      d1,d0      { save quotient.h                       }
  857.    swap        d0
  858.    move.w      d3,d1      { (d0.h = remainder of prev divu)       }
  859.    divu        d2,d1      { divide dvd.l by dvs                   }
  860.    move.w      d1,d0      { save quotient.l                       }
  861.    clr.w       d1         { get remainder                         }
  862.    swap        d1
  863.    bra         @LDiv7     { and return                            }
  864. {;== case 3) divisor > 16 bits (corollary is dividend > 16 bits, see case 1)}
  865. @LDiv3:
  866.    moveq.l     #31,d3     { loop count                            }
  867. @LDiv5:
  868.    add.l       d1,d1      { shift divident ...                    }
  869.    addx.l      d0,d0      {  ... into d0                          }
  870.    cmp.l       d2,d0      { compare with divisor                  }
  871.    bcs         @LDiv6
  872.    sub.l       d2,d0      { big enough, subtract                  }
  873.    addq.w      #1,d1      { and note bit into result              }
  874. @LDiv6:
  875.    dbra        d3,@LDiv5
  876.    exg         d0,d1      { put quotient and remainder in their registers}
  877. @LDiv7:
  878.    tst.l       d5         { must the remainder be corrected ?     }
  879.    bpl         @LDiv8
  880.    neg.l       d1         { yes, apply sign                       }
  881. { the following line would be correct if modulus is defined as in algebra}
  882. {;   add.l   sp@(6),d1   ; algebraic correction: modulus can only be >= 0}
  883. @LDiv8:
  884.    tst.w       (sp)+      { result should be negative ?          }
  885.    bpl         @LDiv9
  886.    neg.l       d0         { yes, negate it                      }
  887. @LDiv9:
  888.    move.l      a1,d3
  889.    move.l      a0,d2
  890.    move.l      (sp)+,d5
  891.    move.l      (sp)+,d4
  892.    rts                    { en exit : remainder = d1, quotient = d0 }
  893. @zerodiv:
  894.    move.l      a1,d3      { restore stack...                        }
  895.    move.l      a0,d2
  896.    move.w      (sp)+,d1   { remove sign word                        }
  897.    move.l      (sp)+,d5
  898.    move.l      (sp)+,d4
  899. @zerodiv2:
  900.    move.l      #200,d0
  901.    jsr         HALT_ERROR { RunError(200)                          }
  902.    rts                    { this should never occur...             }
  903. end;
  904.  
  905.  
  906. Procedure LongMod;Assembler;
  907. { see longdiv for info on calling convention }
  908. asm
  909. XDEF LONGMOD
  910.    jsr     LONGDIV
  911.    move.l  d1,d0      { return the remainder in d0 }
  912.    rts
  913. end;
  914.  
  915. {
  916.   $Log: lowmath.inc,v $
  917.   Revision 1.3  1998/07/01 14:28:32  carl
  918.     * LONGDIV bugfixed with signed and modulo
  919.     * LONGDIV bugfix when divisor is less then 16 bits
  920.  
  921.   Revision 1.1.1.1  1998/03/25 11:18:44  root
  922.   * Restored version
  923.  
  924.   Revision 1.5  1998/01/26 12:01:47  michael
  925.   + Added log at the end
  926.  
  927.  
  928.   
  929.   Working file: rtl/m68k/lowmath.inc
  930.   description:
  931.   ----------------------------
  932.   revision 1.4
  933.   date: 1998/01/05 00:33:30;  author: carl;  state: Exp;  lines: +144 -156
  934.   * Several minor bugfixes
  935.   ----------------------------
  936.   revision 1.3
  937.   date: 1997/12/14 19:03:36;  author: carl;  state: Exp;  lines: +11 -12
  938.   + uses now correct register conventions
  939.   ----------------------------
  940.   revision 1.2
  941.   date: 1997/12/01 12:37:20;  author: michael;  state: Exp;  lines: +14 -0
  942.   + added copyright reference in header.
  943.   ----------------------------
  944.   revision 1.1
  945.   date: 1997/11/27 13:55:47;  author: carl;  state: Exp;
  946.   Internal RTL math functions for m68k.
  947.   =============================================================================
  948. }
  949.